function crank_step

%  Solves the heat equation for various M values
%       diff(u,x,x) = diff(u,t) + f(x,t)   for xL < x < xr, 0 < t < tmax
%  where
%      u = 0  at x=xL,xR  and  u = step at 0.25 & 0.75 at t = 0

% clear all previous variables and plots
clear *
clf

% get(gcf)
%set(gcf,'Position', [1201 277 560 725]);
plotsize(560,725)

% set parameters
N=30;
M=5;
tmax=0.1;
xL=0;
xR=1;

% pick time points (by picking the index)
itotal=3;
it(1)=2;
it(2)=(M+1)/2;
it(3)=(M+1);

fprintf('\n    Solution Computed with N = %3.0f and M = %4.0f\n\n',N,M)

% generate the points along the x-axis, x(1)=xL and x(N+2)=xR
x=linspace(xL,xR,N+2);
h=x(2)-x(1);

% calculate numerical solution
for im=1:3

	% generate the points along the t-axis, t(1)=0 and t(M+1)=tmax
	t=linspace(0,tmax,M+1);
	k=t(2)-t(1);
	lamda=k/h^2;
	
	fprintf('\n    Lamda = %5.2e\n\n',lamda)
	if im==1
		ue=crank(x,t,N+2,M+1,h,k,lamda);
		tt(1)=t(it(1)); tt(2)=t(it(2)); tt(3)=t(it(3));
	elseif im==2
		uee=crank(x,t,N+2,M+1,h,k,lamda);
	else im==3
		ueee=crank(x,t,N+2,M+1,h,k,lamda);
	end;
	M=2*M;
end;

% plot results
xx=linspace(xL,xR,100);
for itt=1:itotal

	% plot numerical solutions
	subplot(3,1,4-itt)
	hold on
	plot(x,ue(:,it(itt)),'-sr','MarkerSize',7)
	plot(x,uee(:,2*it(itt)-1),'-db','MarkerSize',7)
	plot(x,ueee(:,4*it(itt)-3),'-','Color',[0.5 0 0.5],'Linewidth',1)

	% define axes used in plot
	xlabel('x-axis','FontSize',14,'FontWeight','bold')
	ylabel('Solution','FontSize',14,'FontWeight','bold')
	
	set(gca,'FontSize',14)
	
	% plot exact solution
	for ir=1:100
		u(ir)=0;
		for ii=1:50
			an=2*(cos(0.25*pi*ii)-cos(0.75*pi*ii))/(pi*ii);
			u(ir)=u(ir)+an*exp(-pi*pi*ii*ii*tt(itt))*sin(pi*ii*xx(ir));
		end;
	end;
	plot(xx,u,'--k')
	
	box on
	say=['Time = ', num2str(tt(itt))]; 
	if itt==1
		yt=0.91;
		axis([0 1 0 1]);
		set(gca,'ytick',[0  0.5 1]);
		legend(' M = 5',' M = 10',' M = 20',' Exact','Location','South');
		%legend(' M = 5',' M = 20',' Exact','Location','South');
		set(findobj(gcf,'tag','legend'),'FontSize',12,'FontWeight','bold'); 
	elseif itt==2
		yt=0.91;
		axis([0 1 0 1]);
		set(gca,'ytick',[0  0.5 1]);
	else
		yt=0.45;
		axis([0 1 0 0.5]);
		set(gca,'ytick',[0 0.25 0.5]);
	end
	text(0.81,yt,say,'FontSize',14,'FontWeight','bold')
	hold off
end;
say=['Heat Equation: exact vs C-N method when u(x,0) has jumps'];
title(say,'FontSize',14,'FontWeight','bold')


% c-n method
function UC=crank(x,t,N,M,h,k,lamda)
UC=zeros(N,M);
for i=1:N
	UC(i,1)=g(x(i));
end;
a=-2*(1+lamda)*ones(1,N); b=lamda*ones(1,N); c=b; z=zeros(1,N);
a(1)=1; c(1)=0; a(N)=1; b(N)=0;
for j=2:M
	for i=2:N-1
		z(i)=-lamda*UC(i+1,j-1)-2*(1-lamda)*UC(i,j-1)-lamda*UC(i-1,j-1)+k*f(x(i),t(j))+k*f(x(i+1),t(j-1));
	end;
	UC(:,j) = tridiag( a, b, c, z )';
end;

% subfunction f(x,t) '
function q=f(x,t)
q=0;

% subfunction g(x)
function q=g(x)
q=0.5*(sign(x-0.25)-sign(x-0.75));

% tridiagonal solver
function y = tridiag( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end;
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end;

% subfunction plotsize
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);